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1.0 SUMMARY 


The current practice for applied analysis in the aerospace industry is to use 
specially selected combinations of coupled zonal models - inviscid, shear 
layer, etc. - to approximate the field equations of fluid mechanics for 

various applications. The zonal models involve systems of ordinary and 

partial differential equations (field equations). These equations are 
simulated with numerical methods which possess two types of numerical errors - 
residual errors and truncation errors. Residual (solution process) errors 
arise due to insufficient iterations of implicit algebraic equations by 
relaxation or the inversion of ill-conditioned matrices. Truncation (grid 
related) errors arise due to the selection of the grid, the gnd-related 
algebraic equations and the associated boundary and initial conditions. 
Residual errors and truncation errors can be very significant. Assessing the 
effects of numerical error on the modeling of the field equations is a tedious 
and expensive process involving parametric cycling through various tolerances 
on residual error and various grid densities and distributions. Grid length- 
scale control to properly resolve shock and shear layer singularities is 

unavailable except for specialized cases. Because of these difficulties, 

numerical error effects are not commonly examined or controlled with precision 
using conventional methods of numerical analysis, and this can lead to a 
misinterpretation of computed results. 

Efficient solution of the equations of fluid mechanics requires the 
availability of adaptive mesh generation and numerical error assessment 
methods to define numerical errors and guide the grid adjustment process. The 
overall goal of this research program is the development of these methods. 

The objectives of the work reported herein were to begin development of 
algorithms to define error norms (for use as resolution monitors) for 

numerical solution of PDE's and to begin development of a multi-level adaptive 
grid technique for application to the solution of the various equation sets 
used to model fluid flow. The present work is an initial exploratory 
investigation of resolution monitors and adaptive grid technology. 
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The approach was as follows. The literature on error assessment and 
multi -grid methods was briefly reviewed. From this, conventional error 
assessment methods were defined and are briefly described. Three variable- 
order-accuracy defect-correction approaches were identified. The 
one-dimensional (1-D) incompressible potential equation was selected as a test 
bed to investigate error assessment and multi-grid methods. A test problem, a 
channel with a constriction, was selected for which analytic solutions were 
available. The equation was solved for the test problem using point 
relaxation for a range of mesh densities and distributions and the various 
error assessment techniques were evaluated. The test problem was also solved 
using point relaxation and a multi-grid scheme and the characteristics of the 
multi -grid method were evaluated. 

One result is that multi-grid schemes are promising as a basis for developing 
resolution monitors and adaptive grid techniques. Brandt's methodology 
appears to be the most suitable approach to adaptive-grid-control. The 
present study suggests that the multi-grid technology is conceptually 
straightforward to apply to conventional computer codes which solve elliptic 
problems. A second is that for the test problem, reliable estimates of the 
maximum glooal error were obtained from solution output for a number of grid 
levels. From the work completed, it is expected that substantial improvements 
are possible for assessing and controlling numerical errors. A third result 
is that significant improvement for efficient residual error control was 
demonstrated with the test problem. Further work is, however, required to 
develop the three key elements: (1) error norms to guide grid adjustment for 
truncation error control, (2) methods for efficient residual error control 
(relaxation schemes that work well with irregular mesh intervals), and (3) 
adaptive mesh structures based on these error norms. The efficient 
interaction of these three key elements is necessary to obtain adaptive 
solution of the Navier-Stokes equations. 

A follow-on research program is recommended which addresses development of the 
three key elements defined in the work reported herein and noted above. The 
development of these elements would occur simultaneously, utilizing a series 
of research computer programs of increasing complexity. 

Work reported herein was supported by NASA contract NASl-16408 and Boeing IR&D 
funds. 
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2.0 INTRODUCTION 


The Navier-Stokes equations with continuity, energy and state equations are 
the accepted analytical model for fluids whose constituative properties are 
Newtonian. They apply to the flight envelope of most aircraft in the earth's 
atmosphere and to all steady and unsteady laminar and turbulent flow processes 
which influence the performance of these aircraft. The development of 
numerical techniques to model the Navier-Stokes equations has been 
revolutionary in recent years and the pace is accelerating. 

The Navier-Stokes equations are always simplified to related but different 
partial differential equations (POE). These simplified PDE systems are chosen 
to model the essential properties of the Navier-Stokes equations for the 
intended application. In the solution of many of these flows, it is difficult 
to sort out modeling errors from numerical errors. An obvious example is the 
use of the Reynolds averaged Navier-Stokes equations for turbulent flow. 

Selecting appropriate PDE systems depends upon having an understanding of the 
PDE solution properties. These are defined through analytical and numerical 
methods. PDE modeling depends upon understanding the numerical results 

including distortions induced by numerical errors. It is important to control 
these numerical distortions to within tolerances that are consistent with the 
intended application. The cost to achieve a given level of accuracy is also 
important. The cost/benefit relation must be readily accessible, otherwise 
unrealistic levels of accuracy may be demanded without real benefits or 
insufficient accuracy may occur with misleading results. 

The present effort is directed at the physics of smooth steady flows with 
interacting singularity regions such as shocks, boundary layers, free shear 
layers, flame fronts and contact surfaces. Existing solution techniques for 
the equations describing these flows are usually unable to control the length 
scales of the mesh in these singularity regions sufficiently to accurately 
resolve these flow features because they are inefficient at the necessary grid 
scales. As a consequence, computed results are of low accuracy. 
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Efficient numerical modeling of these equations with systems of algebraic 
equations for a grid is difficult because of residual errors in the solution 
process and truncation errors. Residual errors occur because of the Gibbs' 
effect (wiggles or high frequency oscillations in the solution) and the 
stiffness of the equations (acoustic, diffusive and convective stiffness). 
Stiffness is defined as slow convergence toward the target of zero residual 
errors. Truncation errors are due to the grid selected, the grid related 
algebraic equations solved, and the boundary and initial conditions - all used 
to approximate the field equations in an analysis domain. These problems are 
compounded by the grid requirements of singular flows; grid length scales are 
needed near singularities that vary by orders-of-magnitude from those required 
in regions of low gradients in flow properties. 

Algorithms for numerical solution of the Navier-Stokes equations are sought 
which address simultaneously the requirements for 

a. grid related algebraic solution procedures for improved reduction of 
residual errors 

b. error monitors that efficiently assist the grid adjustment process 
and optimize the residuals relative to the truncation errors 

c. solution procedures which are less sensitive to mesh and permit grid 
nesting. 

Adaptive grid control is the interaction of these elements to reduce numerical 
error. Manual and automatic processes can be used to implement adaptivity. A 
balance between computer code development time and computer code user 
complexity must be kept in mind. 

The multi-level adaptive grid procedure^^^ produces truncation error 

estimates as a by-product of the solution procedure which could be used as a 
resolution monitor. This procedure is thus especially attractive as an 
approach to the development of automatic PDE solvers which control numerical 
errors to a prescribed tolerance. 

Work reported herein was supported by NASA Contract NASl-16408 and Boeing IR&O. 
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2.1 


THE OBJECTIVES OF THE STUDY 


The first goal of the research is the development of numerical error 
assessment methods for use as grid resolution monitors. The second goal is 
adaptive mesh generation methods to refine the grid locally where indicated by 
the resolution monitor. The third goal is to improve the efficiency of 
residual error control with non-uniform meshes. It is expected that the 
availability of this technology will provide significant improvement in the 
numerical solution of the PDE's of fluid mechanics. Substantial work is 
necessary before this will be possible. The objectives of the present work 
were to begin development of algorithms to define error norms and to begin 
development of multi-level adaptive grid techniques. The present work is 
an initial exploratory investigation of resolution monitors, grid adjustment 
methods, and residual error control efficiency. 

2.2 THE TECHNICAL APPROACH 

Two overall strategies are being used to guide the development of resolution 
monitors and adaptive grid methods. The first strategy is to define as early 
as possible all of the elements necessary to the development of the desired 
technology and to address these simultaneously. The second is to use simple 
one-dimensional numerical "test beds" to define and develop the necessary 
technology elements. Once the technology elements are defined and developed, 
extension of the technology to test bed codes for 2-D and 3-D PDE's of fluid 
mechanics should be relatively straightforward. With this background, the 
error-norm adaptive grid technology can then be applied to codes for efficient 
solution of fluid flow analysis problems. 

Specifically for the work reported herein the detailed technical approach was 
as follows: 

1) The literature on error assessment and multi -grid methods was briefly 
reviewed. Error sources were identified. Section 3.1, and a brief 
mathematical description of these is presented in Sections 3.1.1 and 
3.1.2. Control of numerical error with filtering and damping and the 
necessary interaction with error assessment methods are described in 
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Section 3.1.3. The problem of efficient residual error control is 
discussed in Section 3.1.4. 

Conventional approaches to error assessment and control are discussed in 
Section 3.2. Three approaches were identified; the conventional 

certification process. Section 3.2.1, the engineering approaches. Section 
3.2.2, and the error norm approaches. Section 3.2.3. 

Four error norm approaches to numerical error assessment were identified. 
Conventional error norms are defined in Section 3. 2. 3.1. A Taylor series 
error monitor approach is described in Section 3. 2. 3. 2. Variable order 
accuracy algorithms for error assessment are described in Section 
3. 2. 3. 3. Multi-grid error norms are then described in Section 3. 2. 3. 4. 

2) Solution of the one-dimensional potential equation was selected as a test 
bed to investigate error assessment and multi -grid methods. Section 4.0. 
Numerical solution of the potential equation using point relaxation is 
described in Section 4.1 and using point relaxation and a multi -grid 
procedure, in Section 4.2. 

3) The test problem was solved using the point relaxation and the multi -grid 
scheme as described in Section 4.3. 

4) The various error norms were evaluated as described in Section 4.4, and an 
adaptive mesh example is presented in Section 4.5. 

5) Computed results were evaluated and are discussed in Section 5.0. 
Conclusions and recommendations for further work are suggested in Section 
6.0 and 7.0. 
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3.0 NUMERICAL CONSIDERATIONS 


3.1 TYPES OF ERROR SOURCES 

The numerical methods for solving PDE's have two types of error sources: 
grid-related solution process (residual) errors and grid-placement related 
(truncation) errors. Explicit marching techniques (time or space) do not have 
residual errors unless there is an implicit equation imbedded in the marching 
scheme. Residual errors occur when implicit equations are solved by explicit 
marching techniques, specialized relaxation schemes or matrix inversion 
processes. Only ideal difference schemes have no truncation-error effects. 
Practical flow analysis tools are not ideal. Truncation-error and 

residual -error effects produce many curious phenomena which must be understood 
in order to develop reliable error norms. The error norms must be able to 
detect any spurious or peculiar numerical phenomena. 

3.1.1 Matheaatical Description 


Let 

LU = 0 3. 1.1-1 

represent the PDE system of interest. 

In discretized operator notation, equation 3. 1.1-1 is 

L^U^ = T^-R^ 3. 1.1-2 

where L^ is the discretization operator, U^ is related to the discretized 
dependent variable vector, and T^ is the local truncation error and R^ is 
the local residual error for each cell of the analysis domain. The grid 
structure index, I, is related to choices of the grid density distributions 
for each selected trial grid where in general the computational grid features 
coupled conformal grids with grid nesting in sub-regions. Ideally the grid 
adjustments are made in some pattern that tend toward a limiting grid 
configuration or 'goal grid.' Thus each unique grid shape is represented by 
-an integer value of I. The 'goal grid' is assigned the index 1^, which is 
known once the numerical error has been constrained to the desired bound. 
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An ideal or perfect difference scheme for 3. 1.1-2 is one in which the local 
truncation error does not contaminate the variables of interest, such as 
velocity, density, pressure, etc. Only the residual errors impact these 
variables. Thus, the user specifies exactly the locations in the geometry at 
which values of these variables are desired. With residual error control 
within adequate bounds, the accuracy of the result is insured within selected 
1 imits. 

Non-ideal difference schemes are defined as those in which the local 
truncation error and residual errors simultaneously influence the value of the 
decoded variables. Except for specialized difference schemes for model 
problems, difference schemes for conventional applied analysis are non-ideal. 
Control of both error sources is addressed herein with emphasis on controlling 
the truncation error. This subject is closely related to the problem of 
proper grid adjustment from an initial state to the 'goal grid' state, with 
proper residual control during the grid-adjustment process. 

3.1.2 Truncation (Grid Related) Errors 

Truncation errors are due to the selection of the grid, the grid related 
algebraic equations, and the boundary and initial conditions which approximate 
the field equations of interest. 

The local truncation error is formally defined as the magnitude that the 
left-hand side of (3. 1.1-2) yields for each cell when the 'goal grid' solution 
is interpolated (restricted) to a trial grid difference equation minus the 
interpolated value of the left-hand side of (3. 1.1-2) from the 'goal grid'. 
It is set at zero in conventional representations of (3. 1.1-2) for all values 
of I. In the defect- and deferred-correction methods discussed in Sections 
3. 2. 3. 3 and 3. 2. 3. 4, truncation error estimates are used to correct for 
truncation error effects. The right-hand side of equation (3. 1.1-2) can also 
be a high-order accurate truncation error correction. 
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This equation exhibits many of the features of the convection terms of 

Navier-Stokes equations but it is simple enough that the powerful methods of 

linear analysis apply. Such analysis indicates that all the numerical methods 

that have been applied to model this equation have the following properties: 

( 2) 

Control of the phase errors, amplitude error and the Gibbs-effects errors 
within selected accuracy bounds depends upon the proper grid density selection 
per period of propagation for wavelengths of interest and with properly 
designed Gibbs-effects filters. Truncation errors have an accumulative effect 
upon the accuracy. 


Gibbs-effects errors are high frequency oscillations near the juncture of 

sharp changes in gradients. For the linear case, the Gibbs' error wavelength 

is about four, seven, sixteen and infinite (zero Gibbs effect) mesh intervals 

for eighth, fourth, second, and first order accuracy convective difference 

schemes, respectively. For nonlinear cases such as near shocks, the wave 

length of the Gibbs-effect error is about two mesh intervals for most schemes 

( 2 ) 

of all orders of accuracy above two. Wavelength smoothing is very 

effective for controlling the amplitude of the Gibbs effect to within chosen 
bounds without introducing global damping. Low order accurate Gibbs-effects 
filters are especially useful for providing the global damping that is 
necessary for removing transient waves from certain types of relaxation 
processes. When properly tuned for the most effective use of the grid, the 
low order accurate filtering devices while serving to adequately damp the 
global waves may not provide sufficiently for the control of the Gibbs' 
oscillations within desired bounds. Wavelength filtering in addition to the 
global damping is well suited to the control of the Gibbs' oscillations within 
desired bounds because the smoothing can be localized as desired. 

Intelligent use of smoothing is one of the central difficulties of modeling 
transient and steady state analysis involving mixed elliptic/hyperbolic 
equations. Where analytical solutions are unavailable, the error norms for 
analysis of required accuracy must account somehow for the phase, amplitude 
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and Gibbs-effects errors. The question whether the solution processes must be 
free of Gibbs-effects errors should be treated in future work? 

Dissipative and non-dissipative convective difference schemes permit expansion 
shocks, artificial gross separation, etc. to form under certain conditions. 
Artificial diffusion is added to eradicate the expansion shocks. Tuning the 
artificial diffusion coefficients for peak accuracy is troublesome. The goal 
of the tuning process is that the artificial diffusion must decay globally and 
locally with mesh refinement so that the accuracy of the sonic line shape and 
position increase with mesh refinement. The error norm used must ensure this. 

3.1.4 Residual Errors 

The numerical modeling of the potential equation or of the Navier-Stokes 
equations leads to algebraic forms in which the propagation of signals is 
retarded as the grid density increases. This stiffness problem leads to 
inefficient reduction of residual errors among the simultaneous system of 
algebraic equations, often leading to exponential or power function decrease 
in convergence as the grid density increases. Typically diffusive stiffness 
is evident in potential flow codes. In Navier-Stokes codes, three stiffness 
problems can appear simultaneously or separately - acoustic, convective, or 
diffusive stiffness and can be aggravated by non-uniform grid. Any or all of 
these factors can undermine the convergence rate severely and can enlarge 
errors substantially. The error monitors must be designed to detect slow 
convergence or inefficient residual error control. For some error norms this 
is a severe requirement. When the local residual error is reduced to some 
fraction of the local truncation error, further reduction of the residual is 
not cost effective; determination of this fraction is a subject for further 
study. 


3.2 ERROR ASSESSMENT AND CONTROL 

Three approaches are used to assess the accuracy of PDE modeling. The first 
is to construct an error difference table with solutions on trial grids of 
various mesh densities and distributions. The second is the use of auxiliary 
information such as experimental data and related analytical solutions for 
various regions of the analysis domain, and the third is the use of the local 
truncation and local residual error estimates associated with suitable error 
norms and error bounds. These three are referred to as the conventional, 
10 



engineering (analogical) and the error-norm approaches to certifying the 
accuracy of PDE modeling. Direct numerical evidence of the accuracy of the 
PDE modeling results from the conventional and error-norm approaches. These 
methods have their origins in numerical analysis technology. In the 
engineering approach no attempt is made at achieving direct evidence. 
Inferential reasoning is predominately used. Discussions of this are given in 
the following sections. 

3.2.1 Conventional Certification Process 

The process of numerical error assessment with conventional PDE modeling 

(3) 

techniques is the following.' ' A solution of finite difference equations 
(simultaneous system of algebraic equations) for a specific discretization of 
the analysis domain is generated for different choices of grid density and 
grid distribution in the analysis domain. It is common to use a sequence of 
grids of the same grid distribution that differ in grid count in each 
independent variable direction by factors of two — 2, 4, 8, 16, 32, etc. The 
coarser grids can be generated by deleting every other point of the finer 
grids. The effects of the choice of grid distribution are examined by 
choosing sequences of grids which have different mesh distributions. The data 
from all of these solutions of the grid-related equations is organized by 
constructing an error difference table. Solution differences are posted in 
order of the coarse-to-fine grids for each grid sequence. The solution 
differences are generated by subtracting the values of adjoining pairs of grid 
solutions of the dependent variables at all physical locations in the analysis 
domain that correspond to the grid coordinates of a grid of a selected 
intermediate density. Interpolation is used to relate other grid solutions to 
these selected grid coordinates. As the grid density increases the 
differences should decay approximately* according to the formal order of 
accuracy for some selected mesh distribution. Iterative adjustment of the 
grid distribution and density is made until this type of error decay is 
realized. If this occurs, extrapolation may be used to solutions at infinite 
grid density and reliable estimates of the maximum global error may result. 


*Error decay according to the formal order of accuracy is expected globally 
but not locally in regions of singularity. 
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The preceding process appears to work best on the modeling of parabolic and 
elliptic equations in smooth domains with smooth boundary conditions. For 
mixed elliptic/hyperbolic systems, erratic results may occur due to unresolved 
singularity regions and/or poor residual error control, and/or Gibbs' error 
effects. 

A key feature of this method is that grid adjustments are made in some pattern 
that tends toward a limiting grid configuration. A way to think about this is 
to define a goal grid to which the selected grid sequences must evolve. The 
'goal grid' is a grid upon which the solution is sought to some specified 
error bound. It should be understood that the 'goal grid' may not be unique 
because of grid initialization, grid generator, and grid-equation solver 
properties. It is assumed that adequate control of the residual error effects 
has been maintained in the process of assessing the truncation error effect. 
This is done by developing a sequence of several solutions on each grid choice 
with various choices of constraints on the residual tolerances that are used 
to terminate the computations for each solution on that grid. Because of the 
need to assess the contamination by residual error, the 'goal grid' may not be 
the grid of greatest density but it will have the correct shape. 

Conventional techniques for developing the data necessary to certify the 
accuracy of numerical modeling procedures are limited by the following factors: 

1) Costs. 

2) Because of (1) above, a very limited number of solutions and thus 
sparse information are usually available from which error estimates can 
be constructed. 

3) Because grid adjustment to control the error within desired bounds is 
cumbersome or impractical, arriving at proper grid configurations in 
mesh density and distribution is also difficult or impractical. 

4) Numerical error during grid refinement may vary erratically, not 
monotonically with grid density. Confusion as to grid adjustment needs 
can result. 
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5) The software is usually not available for conveniently constructing the 
error table. This means that the error assessment process is manpower 
intensive. These factors discourage development of the 'goal grid' 
solution. Without this, the accuracy of the result is unknown; the 
meaning of the result is undefined and useless unless appropriate 
external information (user experience) is applied to the result. 

The certification process with conventional computer codes utilizes error 
information from many grid densities and grid geometries. This is a 
multi-grid process albeit a very cumbersome and inefficient one. For the 
present purposes this technology will be referred to as fixed grid (FG) even 
though it is not, when properly used for numerical error assessment. It is 
called FG because that is the manner in which it is used in applications; the 
error term associated with each number in the output is set to zero and this 
is often ignored in the use of the numbers from the output. 

3.2.2 Engineering Approaches 

Judgement in engineering applications as to what grid should be selected for 
numerically modeling a PDE system is often based upon an exterior body of 
knowledge' ‘ rather than direct numerical evidence. For example, 
boundary-layer analysis can be performed with finite difference and related 
analysis tools. Mesh refinement studies and conventional error norms can be 
used to define the accuracy and grid-choice relationships. Similar studies 
can be performed on free shear layer and inviscid model problems with 
analytical solutions to establish the accuracy and grid-choice relationships. 
The various component features of the physics of the PDE system can be studied 
in this manner. The selection of the trial grid for the PDE system of 
engineering interest can be made upon the basis of the physics that is 
expected in each flow region of the analysis domain. The grid selection in 
the various flow regions of the analysis domain can reflect the desired 
accuracy that is required locally and globally for the purposes of the 
analysis. Interpretation of the results of numerically modelled PDE systems 
involves qualitative aspects of the solution. Inspection is used to insure 
that solution features such as wall shear stress, wall boundary layer, free 
shear layer, inviscid structure, or shock structure characteristics occur 
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where they are expected. The success of this approach depends on the 
knowledge and skill of the analyst and the time allotted for the analysis. 
Important physical processes may be inadvertently ignored because numerical 
errors mask solution properties. For example, artificial diffusion can be 
interpreted as turbulent diffusion. 

Another approach to grid selection is used in engineering applications as well 
and sometimes augments the above approach. It involves the comparisons of 
computed and experimentally measured flow properties. Grid and local and 
global numerical smoothing adjustments are used to generate favorable 
agreement between the computed and measured flow properties. Where 
experimental data is available, this approach is preferred to those that are 
described above and it encourages the use of analysis for predictive purposes 
where favorable agreement occurs. 

In engineering approaches, direct numerical evidence is not used to understand 
the nature of the numerical properties; rather inferential and analogical 
reasoning are used. 

3.2.3 Error Norm Approaches 

Four possible methods are considered here. The first method is a brief review 
of conventional error norms; the second is the use of truncated Taylor series 
expansions for an error monitor; the third is the use of variable order 
accuracy algorithms; and the fourth is the multi -grid approach. These are 
discussed below. 


3. 2.3.1 Conventional Error Noms 

( 4 r T 

standard error norms include: 

El = C z|((!)^ - (t)^^|]/N 
E 2 = Cz((|)^ - 

E = [|((t)^ - 
max I ^ ^ ' 

(p = ampl i tude 



3 . 2 . 3 . 1 - 1 

3 . 2 . 3 . 1 - 2 

3 . 2 . 3 . 1 - 3 
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where (|)^ is the dependent variable of the analytical solution restricted to 
the coordinate location of grid point I. N is the number of grid points, 
(j)^^ is the discretized solution on grid I. The sums denote component sums 
in 2-D and 3-D situations. and are known as the average 

error, the root-mean-square (rms) error and the maximum global error, 
respectively. The norm is best applied to non-singular problems (no 

shocks, shear layers, or geometry discontinuities) with smooth boundary 
conditions although it can be used if the regions of steep gradients are 
excluded. Ej^ and E2 are used for singular perturbation problems including 
the steep gradient regions. The use of these definitions as written depends 
upon knowing (j)^ which is not available for cases of general interest. 
However, modified forms of these relations may be considered. 

For example, equations 3. 2. 3. 1-1 through 3. 2. 3. 1-3 can be given by: 


- (t)^^)l6(Z) ]/N 3. 2. 3. 1-4 

3. 2.3. 1-5 

= CG(Z)1«|)^^ - /(M^\ax 3. 2. 3. 1-6 

where is any of ... is also regarded a 

higher order accurate solution on grid ll. G(Z) is a weighting factor which 
may be used to exclude regions in the analysis domain where the local 
truncation error estimates exceed a selected threshold value. N is the number 
of cells in the summation. These error norms should be of general use for 
understanding the role of residual and truncation errors, especially when used 
with the error norms discussed below. 

There are error estimators in addition to the ones listed above which should 
be considered in any application of interest. For example, certain components 
of PDE systems require that kinetic energy, mean vorticity, total pressure, 
entropy, mass, etc. should be conserved. Discretized forms of these integral 
relations could be constructed which are relevant to the application of 
interest. 
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3. 2. 3. 2 Taylor Series Error Monitor 

The idea of using truncated Taylor series expansions for an error monitor has 
been explored and conclusions about the cost effectiveness of such a 
development are summarized as follows. Consider a simple harmonic function 
such as the y=sin x. From calculus, the derivatives of all orders exist. The 
odd-order derivatives vanish when the modulus of the even-order derivatives 
are maximum. At the least, pairs of odd-order and even-order derivatives 
would have to be approximated by finite difference expressions to prevent 
spurious decay of higher-order terms in any proposed error monitor. The 
second requirement is that, in regions of the grid with long wavelength 
variations of a dependent variables of interest, it must establish how many 
pairs of odd and even orders are required to certify the numerical accuracy of 
the PDE modeling. The third requirement is that, since high frequency 
oscillations (two mesh interval limit) can excite two mesh interval 
oscillation on all discretized derivative of the order of two and above, this 
information must be used somehow to sort Gibbs errors, geometric 
discontinuities, etc. While the Taylor series error monitor idea is 
intriguing, it has conceptual and implementation problems. The next section 
describes a related but perhaps better approach to error assessment. 

3. 2.3.3 Variable Order Accuracy Algorithms for Nianerical Error Assessment* 
Rather than using only mesh refinement to assess numerical error, it is 
conceptually reasonable to define numerical algorithms of variable order of 
accuracy from which direct error estimates emerge for a given grid. A way to 
think about this is to define a sequence of higher order accurate discretized 
solutions on the same grid each of which satisfy a certain smoothness 
property. The higher order accurate solutions on the grid are regarded as the 
trial analytical or reference solutions. The lowest order accurate solution 
is regarded as the approximation. The error estimate is performed with 
conventional error norms. The error norms of Section 3. 2. 3.1 are slightly 
modified for error assessment by using higher order accurate solution data in 
place of the 12 finer grid data. Iteration can be used to determine how high 


* Defect- and deferred-correction schemes are closely related because error 
estimates are used to improve the efficiency of residual error control. 
Both methods are multi grid methods. 
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the order of accuracy must go for reliable error estimates. Grid adjustments 
are needed to provide the highest order accurate scheme with sufficient grid 
to meet desired accuracy bounds. 

To implement the above process, a convergent variable order of accuracy 
algorithm approximating the PDE system of interest must be available. For the 
boundary layer equations, Wornom^^^ presents an approach. Forester^^^ 
suggests a method for mixed hyperbolic/parabolic PDE systems using odd-degree 
splines although other basis functions can be used. The order of accuracy of 
this scheme is selected simply by choosing the appropriate coefficient matrix 
that is associated with the desired accuracy. With this method, orders of 
accuracy six, ten, fourteen, nineteen, etc. are possible and they are 
activated by an input control function. Interpolation is used to initialize 
successive solution processes of higher accuracy. 

Lindberg^®^, Pereyra^^^ and Stetter^^^ have developed an approach which 
is related to the one described above but in which certain simplifications are 
implemented. Basically Lindberg uses a low order accurate approximation to 
the PDE system with a non-zero right-hand-side term that is lagged in the 
iteration process. This perturbation operator is a high-order-accurate 
discretization (usually fourth order) of the PDE system. The data for the 
perturbation operator is derived from the low order of accuracy discrete 
equation system solution. The perturbation equation is repeatedly solved 
until some convergence criterion is satisfied. The output of the perturbation 
equation is a correction parameter. It is added to the low order of accuracy 
solution to yield a fourth order accurate solution at convergence. Error 
estimates are produced with conventional error norms by using the high-order- 
accurate solution as an exact representation of the PDE system. The low order 
accurate result represents the approximate solution. Lindberg presents many 
examples of the efficiency of this defect-correction procedure to reduce 
residual and truncation errors. Success is not universal with hyperbolic 
problems, however, when some smoothness property is violated. It may be 
inferred from Lindberg' s data that if sufficient smoothing is applied to 
regions of rapid change in gradient, convergent results are achieved. To make 
Lindberg' s scheme useful for hyperbolic systems, empirical studies would have 
to be performed to develop criteria for how local the smoothing must be 
constrained so that the global errors are definable. See Section 5.3 for 
related discussion. 
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The approaches of references (1-8) presume that some sort of mesh refinement 
process toward the 'goal grid' is used. The accuracy of these error 
estimators improves dramatically as the accuracy moves toward the one percent 
error range. They are useful guides for grid refinement when errors are 
larger than ten percent but are not precise in the predicted magnitude. These 
approaches are promising for modeling mixed parabolic/elliptical /hyperbolic 
systems of equations but are difficult to implement when grid nesting is 
required. 

(9) 

Another approach to defect-correction has been proposed by Brandt . This 
approach suggests variable order of accuracy, nested grids, ennancement of 
residual control efficiency and the efficient generation of an error term that 
may be useful for numerical error assessment. This approach has been called 
the multi grid method but it is misnamed since all numerical assessment and 
control schemes are multi grid.* Many variations of this methodology are 
possible. Most examples of this approach^^’®’^^^ are for residual control 
efficiencies near theoretical limits with uniform grid intervals. Irregular 
grids (nested grids) have yet to be widely addressed in this method. 

3. 2. 3. 4 Multi -Grid Error Noras 

Modifications are suggested in Section 3. 2. 3.1 to conventional error norms for 
assessing numerical error schemes. This is applicable to the FAS-MG (Full 
Approximation Storage-Multi -Grid) scheme with one exception. The solution 
data at each grid level must be saved before finer grid levels are invoked in 
the multi-level solution process. This approach insures that the targets of 
residual and truncation errors are zero in the output for the various levels 
of grid. If this is not done, solution errors may appear on coarser grids 
when in fact none exist. 

The FAS-MG scheme estimates local truncation error by interacting coarse and 
fine grid solution data. Brandt suggests that the weighted integral 

E-,. = JG(z)jT(z)|dz 3. 2. 3. 4-1 

* Brandt's approach could be called defect-correction multi-grid or Tau multi- 
grid, TMG. 
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is a measure of global error level, Ej, where G(z) is either zero or unity. 
G(z) represents a weighting to exclude singularities in the analysis domain. 
T(z) is the local truncation error. The discussion in paragraph 3.1.4 
suggests that additional error norms should be considered for the control of 
residual errors. For example, relationships such as 

Ej^ = jlR(z)ldz 
^RT = IlR(z)l/lT(2)ldz 
may be suggested. 

These error norms can be readily generalized for more space dimensions. R(z) 
and T(z) are normalized by some suitable reference quantity so that they 
cannot exceed unity during the solution process. Singularities in the POE 
modeling are detected by T(z) of the order of unity. 


3. 2. 3. 4- 2 

3. 2. 3. 4- 3 
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4.0 KUMERICAL ANALYSIS AND RESULTS 


With reference to the technical approach. Section 2.2, solution of a 
one-dimensional potential equation was selected as a test bed to evaluate 
error assessment and multi-grid methods. Solution of the one-dimensional 

potential equation was selected because of its simplicity, only a single 
dependent variable. The equation was solved using a point relaxation scheme 
alone and then a multi -gridded point relaxation scheme. 

A test problem, a channel flow with a constriction, was selected because 
analytical solutions are available for comparisons with the numerical 
solutions. The test flow was solved for a range of mesh densities and 

distributions and the more promising error assessment techniques were 
evaluated. The test problem was also solved with a multi-gridded point 
relaxation scheme to understand and evaluate the characteristics of the 

multi-grid scheme. A semi-self-adaptive mesh scheme was briefly investigated. 

4.1 SOLUTION OF THE 1-D POTENTIAL EQUATION 

The discretized 3-D full -potential equation for steady incompressible flow is 
restricted to a 1-D analysis tool by deleting the K and L indicies in the Ref. 
10 formulation. The total velocity can be computed in many ways. One 

formulation^^®^ yields values of the total velocity in which the truncation 
errors in the velocity potential do not contaminate the computation of the 
total velocity. 

The forms of the 1-D continuity equation for the purposes of the present study 
are in terms of the primitive velocity component, W, and the velocity 


potential, (j>. They are given by 


(WA)^ = 0 

4.1-1 

= 0 

4.1-2 

U = (Ij 

4.1-3 
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where A is the channel cross section. 

The transformed discretized forms of these equation are given by 


C((W) (D))„^,^- 


‘m > 


4.1- 4 

4.1- 5 

4.1- 6 


= grid level index, I = 1 coarsest grid, two cells in domain range of 
zero and unity for Z/L 



= (2)(D„,^)V(VL„.VL,^,) 



' ^ " ^M+% 



= 1 

Z less than zero 


= Y(z) 

0 < Z < L cubic functions 

V% 

= 1 

Z greater than 1 


VL, 




M 

M 

Ml 


= residual error due to inexact solution process 
= the grid index of cell centers 
= 1, dummy cell for boundary condition data 
= 2, first cell in analysis domain 
= + 1, last cell in analysis domain 
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M2 


= + 2, dummy cell for boundary condition data 

InaX 

<1>M2 " °M1 boundary value 

= (|)2 - We (Vl^ + VL2)/(D2_,^)/2 boundary value 
= grid index of cell face of cell M, 

Mjjjax " index of cell face on cell M 
We = velocity at entrance of the 1-D channel 


Given We, for 2'^^1, and Rj^ for 2£M_^1, equation 4.1-4 is 

solved exactly for Wj^^,^by marching from M=2 to Ml. The analytical solution 
to equation 4.1-1 is recovered by equation 4.1-4 at any selected z coordinates 
if Rj^ is identically zero during the marching process (round-off error 
effects are ignored for practical purposes.) With and Rj^ 

selected equal to zero, equation 4.1-4 is regarded as a perfect or ideal 
difference scheme. Equation 4.1-5 is regarded a perfect difference scheme in 
terms of Wj^ since equation 4.1-4 results from combining equations 4.1-5 and 
4.1-6. 


Relaxation methods can be used to approximate equation 4.1-5. Two relaxation 
equations are 


a * 0 , Richardson scheme 4. 1-7 a 

a = + , Liebmann scheme 4.1-7b 
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= updated 
i[)j^ = old value 



Note that is set to zero when deriving equaions 4.1-7 by the manipulation 
of 4.1-5. 

The FG method of solving equation 4.1-2 is defined by iterating equation 
4.1-7a or 4.1-7b until some test on the residual show that this process should 
be terminated. 

4.2 SOLUTION OF THE 1-0 POTENTIAL EQUATION WITH MULTI-GRID 

The FAS-MG scheme is more elaborate. It is defined as follows. Equations 
4.1-7a and 4.1-7b are modified by inserting a term, RSj^, into the numerator 
so that they read 


a = 0 , Richardson scheme 4.2-la 

a = +, Liebmann scheme 4.2-lb 



RSj^is the local truncation error estimated on grid I relative to grid Ig 
on which RSj^^ is set to zero. l|g is the restriction operator for 
the conversion of data on grid Ig to a usable form on grid I. The subscripts 
in equation 4.2-lc refer to grid I. The process of generating data 
that is necessary to use equation 4.2-lc is provided in an excellent form by 
Brandt (Ref. 8). The variant of this process adapted here (see diagram 1) is 
to start on the coarsest grid (Level 1) with equation 4.1-7b. When the mean 
residual tolerance of 10”^ or lower (see below) is satisfied, standard 
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linear interpolation of the velocity potential to the next finer grid level is 

implemented. Equation 4.1-75 is iterated on this grid (Level 2) until a 

minimum number of iterations (1, 2, 4, 8, 16, 32 were used) are completed. 

Further iterations are required if the ratio of the new and old residuals, the 

eigenvalue, is less than some amount - .82 was used based upon trial -and-error 

experience. The rate of residual error reduction is slowing if .82 is 

exceeded. When this occurs, it is time to switch to the next coarser grid 

(Level 1). This means that the are interpolated from the Level 2 grid 

to the Level 1 grid and Equation 4.2-lb is introduced. Rsj^ is estimated 

from Equation 4.2-lc with Ig replaced by 2. A converged solution to a low 

residual tolerance is produced on Level 1. Velocity potential data for Level 

2 again must be interpolated from existing data at Levels 1 and 2. For this 

purpose, Brandt’s prolongation* operator is used. The details of this 

equation are provided later in this discussion. It may take several Level 1 

and Level 2 cycles before a low residual tolerance is achieved on Level 2. 

Once convergence on Level 2 is achieved, prolongation to Level 3 is used with 

standard linear interpolation of the velocity potential on Level 2. Equation 

4.1-7b is iterated on Level 3 until the minimum number of iterations are 

performed or until the critical eigenvalue is exceeded. Restriction to Level 

2 

2 provides an estimate of RSj^ from Level 3 using the Equation 4.2-lc with 
Ig replaced by 3. Level 2 is iterated with this value until the solution is 
stalled. Restriction to Level 1 provides a corrected estimate of RS^^ 
Cycling between Levels 1, 2, and 3 continues until Level 3 residual tolerances 
are satisfied. Level 4 is then used. This process can go on indefinitely, 
but it is usually terminated at Level 5. See Diagram 1 for a flow chart 
description of the FAS-MG scheme that is used. 

The R[^ term is composed of two components — the residual and the local 
relative truncation estimate. At the finest grid yet arrived at during the 
multi-level solution process, these elements are of equal and of opposite sign 
so that the target of RSjJ is exactly zero. In regions of the coarse grid 
solutions where the truncation error estimate should be zero, it is not zero 
but it has a magnitude proportional to the residual error on the finest grid 
yet arrived at. Its actual magnitude depends on how many grid levels it is 


*Prolongation is defined as interpolation to a finer grid. 
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Brandt’s FAS— MG Scheme (Modified with Error Control on Residual) 
Coarest-to-Finest Grid Approach to L^“0 



DIAGRAM 1 
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removed from the finest grid yet arrived at. A rule-of-thumb is that the 
residual error of the finest grid yet arrived at is about doubled each time 
the grid intervals are doubled. 

For the present studies, two approaches to the selection of residual 
tolerances are examined. The residual tolerance can be chosen either for the 
coarsest grid or the finest grid. In the former method, the coarsest grid 
residual tolerance is reduced by a factor between two to four each time the 
grid size is doubled. Alternatively, if the finest grid residual tolerance is 
chosen, it can be applied to all grid levels. Both methods have been used. 
The latter method is recommended for simplicity in the use of the modified 
error assessment formulas of paragraph 3.2.3 which are discussed in Section 
4.4. 


4.3 THE TEST PROBLEM 

The 1-D test problem involves an analytical geometry of a straight channel 
with a cubic function for a constriction that reverts either abruptly 
step-wise or smoothly to a straight channel. Figure 1 shows the channel 
section shape distribution with respect to the flow direction. Figure 2 snov.s 
the analytical solution restricted to 65 grid coordinates (64 cells) with the 
grid intervals constant. Eleven trial fine-grid sets were used to examine the 
1-D potential solution properties for a) grid with uniform mesh intervals, 
b)grid with uniform mesh intervals in the region of cross sectional area 
variation but with a stretch factor of two in the straight sections, and c) 
grids with uniform mesh intervals in the straight sections but with a stretch 
factor of .80, .85, .90, .95, 1.0, 1.05, 1.1, 1.15, and 1.2 in the constricted 
region where the finest grid is near the abrupt enlargement of the channel 
cross sectional area for stretch factors less than unity. The total number of 
grid intervals for each set of trial grids are 4, 8, 16, 32, and 64, where the 
number of grid intervals in the constriction region are respectively 2, 4, 8, 
16, and 32. FG and MG methods have been applied to generate solutions, shown 
in Figure 2, by the solid line for the finest grid. Also shown in Figure 2 is 
the FG and MG solutions with very high residual error tolerances. Using point 
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relaxation* and sweeping the grid in the flow direction, MG yields a maximum 
global error** of less than 4% in the equivalent*** of twenty-five sweeps of 
the 64 node grid, whereas FG requires over one thousand sweeps of the 64 node 
grid to achieve the same accuracy. The maximum error occurs at the geometric 
discontinuity. Increasing the accuracy by an order of magnitude requires less 
than a factor of three increase in the work for the MG and the FG.' The 
process of solving the problem to greater accuracy can be continued until the 
maximum global error satisfies desired constraints up to round-off error 
effects. The boundary conditions are imposed both on the FG and MG as set 
mass rates of equal magnitude at the entrance and exit cross sections. 

Figure 4 shows the truncation error spectrum for the peak values of the local 
truncation error asymptotically approach nearly the same values including 
T-extrapolation, on the next-to-the-f i nest grid solution. The magnitude of 
these terms are substantial near the discontinuity and, because they form the 
right-hand side of the cell-wise flux balance equations, induce large errors 
in the total velocity profiles that are shown in Figure 3. The coarse-to-fine 
grid correction equation of Brandt very effectively interpolates the Poisson 
type solutions on coarser grids so that the coarser grid solutions mimic the 
finer grid solutions. Standard interpolation (prolongation), 



cannot account on the next finer grid, I+l, for the fact that the right-hand 
side term is significant in the coarser grid solutions. For this reason, 
standard interpolation is not useful and must be replaced by a more elaborate 
interpolation. Brandt recommends (for prolongation) 

“new I '“new I+l “old' “old 


* Both Richardson and Liebmann point relaxation were used. As expected, the 
efficiency of FG and MG are improved with Liebmann relaxation. 
Conclusions about the asymptotic efficiency with mesh refinement holds 
irrespective of the form of the relaxation used here. All of the 
displayed results are with Liebmann relaxation. 

** The global error is defined by equation 3. 2. 3. 1-3. 

*** See the discussion in Section 5.2 for a definition of equivalent sweeps. 
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where ij , is the fine-to-coarse grid interpolation operator and 
I+l 

Ij is the coarse- to-fine grid interpolation operator. This expression 

functions well as illustrated in Figure 5. Linear interpolation is used for 

these operators with weightings of 1/4 and 3/4 for and weightings of 

I ^ 

1/2 and 1/2 for These weightings are derived directly from the 

geometric relationship between the coordinates of the cell centers of the two 
adjacent grid levels whose cell faces coincide at every other cell face. No 
modification of this weighting is used for stretched grid cases whose 
principle effect is to retard the convergence rate by up to one- third for 
cases with stretch factors of .80 and 1.2. 

4.4 ERROR NORM EVALUATION 

The utility of the error norms of paragraphs 3. 2. 3.1 and 3. 2. 3. 4 is examined 
in this section. Errors in the computed velocity are studied with the use of 
the maximum global error estimator and average and maximum truncation error 
estimators. It is expected that similar conclusions would be reached if other 
error estimators of paragraph 3. 2. 3.1 were used. 

The unmodified (analytical reference) and the modified (finer grid reference) 

II 12 

maximum global error norms (^niax’ ^ax ^ Section 3. 2. 3.1 and the 
error norms of Section 3. 2. 3. 4 have been applied to a number of cases of 1-D 
incompressible channel flow that have smooth and abrupt cross sectional area 
changes. A new error norm is defined and is used as well. The results for 
these error norms are summarized as follows. 

II 12 

To illustrate the properties of and E^iaxl ’ output station 

is chosen for which the grid levels — 2, 4, 8, 16, and 32 cells in the 

transition region of the channel geometry. In all cases, the location is 

selected nearest the minimum channel cross section where the largest errors in 

velocity reside. increases in size as the grid size and/or the 

residual tolerances grow. However, this nice behavior does not occur with 
II 12 II 12 

^max * ^ax unique since the output from the various 

grid levels, 12, is used for the reference solution. 
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Let 12 equal grid levels 2, 3, 4, and 5 to examine the maximum estimated 

1112 

global error for the coarsest grid, 11=1. Emtl values increase in 

IliQ A 

size as the reference grid level and/or the finest grid residual tolerances 
II 12 

grow. values are a measure of relative error between solutions 

at different grid levels and as such may have a different sign and level than 

II 12 

Emax* Furthermore, values define the error in the reference 

grid solution, 12, rather than in the approximate solution, II, under 

examination. This conclusion is predicated upon having the residual tolerance 

for the approximate solution within an error bound that is about the same 

magnitude as the reference grid which causes a more accurate coarse grid 

solution than that of the finer reference grid solutions. This result holds 

strictly only for a perfect difference scheme. For nonperfect difference 

schemes, local truncation error will likely dominate the coarse grid 

solutions. However, it is possible that peculiar local truncation error may 

produce smaller real error in coarse grid solutions. Therefore, it is 

necessary to use additional information to determine if the error indicator 
II 12 

■'s a measure of coarser grid error. One method of determining 
max j. J 2 

this is to examine the behavior of as the residual tolerance is 

max 

reduced. For a nonperfect difference scheme, E„,’ should reach a 

max 

fixed value for some range of residual tolerance. If so, the E^^j^ 
indicator is measuring the coarser grid error. Otherwise, residual error 
effects are dominant. 


:3,5 


The behavior of is illustrated for E^’^, and 

« M 0/1 illuA iiiQA f llldA 

predicts the error in the level 4 solution to about 
thirty percent accuracy of the true error. predicts the error in 


the level 5 solution to about a five percent accuracy of the true error on 


level 5. 


^ax 


predicts the error in the level 5 solution to about a 


one-half of one percent accuracy of the true error. These statistics are for 
a residual error of 10"® at all grid levels. This residual tolerance 
reflects about a one tenth of one percent true error range for the level 5 
solution. For the ten percent true error range of level 5 solutions, the 
precision is reliable to better than an order of magnitude. This provides a 
useful guide for grid density adjustments. 
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II 12 

The behavior of E_,’ is spurious when the residual is controlled on 

each grid level to yield the true error of the same magnitude on each grid 

level. Brandt recommends residual error control in this fashion which renders 
II 12 

useless for ideal difference schemes. For this reason and 

Ilia A 

because nonideal difference schemes may locally, in certain cases, behave like 
an ideal difference scheme, it is suggested that it is better to control the 
residual to the same level on each grid level even though this guideline may 
not yield peak computatonal efficiency exclusive of the error assessment costs. 

Application of Equation 3. 2. 3. 4-1 with G(z) set equal to unity yields the 
average value of the local truncation error. The average and maximum values 
of the local truncation error are examined for utility in error assessment. 
When scaled by the channel cross section, these quantities are converted into 
average and maximum velocity perturbations, respectively. The scale factors 
are the average and the minimum channel cross sectional areas, respectively. 
These velocity perturbations have been correlated to the maximum global error 
in the fluid velocity for various grid levels that contain the nonzero R.H.S. 
on all but the finest grid. The average velocity perturbation predicts 
conservatively the right order of magnitude that the MG generated maximum 
velocity must be corrected by in order to estimate the error in the MG 
output. This result applies to smooth and discontinuous channel shapes with 
uniform grid intervals. The use of the maximum value of the local truncation 
error appears to be best reserved for locating discontinuities in the solution 
since it tends to overestimate the solution error levels in all grid levels. 
This error estimator will be especially useful to identify regions in the 
analysis where the mesh interval is too large irrespective of how coarse the 
grid is in the neighborhood of discontinuity in dependent variables such as 
temperature, pressure, tangential velocity, etc. For the immediate future 
this error estirator has application to all existing flow codes. The 
installation of the FAS-MG scheme in these codes would be useful just for that 
purpose alone. Considerable cost saving could be realized by using this 
method of guiding the grid adjustments. 

Only one exact method of evaluating the maximum and local global error has 
been found for the perfect difference scheme. The sums of the same-sign local 
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truncation error are converted into velocity perturbations by the use of 
Equation 4.1-4. The sign content of these terms oscillate to produce 
non-smooth corrections. Starting at the edge of the analysis domain, these 
sums are added one-by-one. This velocity is called the modified exact 
velocity. At each point where the sign change occurs in the velocity 
perturbation, the error between the grid solution and the modified exact 
velocity is computed and saved. All of these errors are sorted until the 
largest is found. The largest value is the maximum global error on the finest 
grid level to within round-off error effects which means that it has extreme 
accuracy. Roughly the error is bounded by one part in ten to the ten on the 
CDC CYBER 175 computer. 

4.5 ADAPTIVE GRID EXAMPLE 

In this section, uses of the local truncation error estimates for grid 
adjustment are discussed. A simple example of semi -adaptive grid refinement 
is shown in Figure 6 in which grid compression toward the region of high local 
truncation error is used. Iterative grid compression is continued until a 
condition of the maximum normalized local truncation error is less than .08. 
Semi-adaptive grid compression is implemented in the interval 0 Z/L £ 1 
by iteratively decreasing the grid stretch factor from an initial value of 1.2 
in steps of .05. As expected, the tolerance on the maximum local truncation 
error is not satisifed as long as an exact step-wise discontinuity is enforced 
at a Z/L equal to unity. With a cubic transition function in the interval 
31/32 ^ Z/L £ 32/32 which has a slope continuity with the remaining 
channel geomtry, local truncation error reduction results with grid 
refinement. Figure 6 shows the results of the analytical solution and 

solution with a grid contracted toward Z/L equal to unity. Over an order of 

magnitude reduction in the local truncation is readily achieved with a 

contraction ratio of .85. The lower the magnitude selected for the stopping 
criteria, the more grid is compressed into the region of the abrupt geomtry 
change. Eventually this approach starves the remaining domain of the analysis 
of sufficient mesh to satisfy the selected maximum local truncation error 
tolerance. Therefore a preferred strategy involves sub-dividing the region of 
small length scale, 31/32 < Z/L < 32/32, with a uniform grid of varying 
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number of grid points. It is easy to implement. It is regarded also as 
semi -adaptive. A 'fully' adaptive strategy requires labeling each cell of a 
trial grid with a special flag that designates cells with a local truncation 
error that exceeds a selected threshold value. Cells so flagged may be 
sub-divided by nesting compressed grids or by uniform interval grid 
embedding. It is expected that the rapid grid-interval changes may produce a 
growth in local truncation error in that region. If this occurs, criteria 
must be developed for the control of the rate of the grid interval varations 
or the meaning of the truncation error reassessed. 'Fully' adaptive MG 

strategy only requires that iterative work to reduce the truncation error be 
applied to the flagged cells. This approach may be more efficient, 'fully' 
adaptive and more computer programming intensive than the semi -adaptive 
strategies. This approach appears to be practical to program for machine 
computations for multi-dimensional numerical analysis. 
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5.0 DISCUSSION 


5.1 ERROR ASSESSMENT 

For the test problem, local normalized truncation error estimates of the order 
of unity seem to indicate the region in which grid adjustment (mesh density or 
distribution) should occur or the region in which the geometric representation 
of the boundary of the analysis domain may need modification. The local 
truncation error estimates in themselves cannot distinguish the cause of large 
local error or whether the results of the analysis are adversely affected. 
Therefore additional information must be associated with the local truncation 
error estimates to make them useful. The behavior of the solution in high 
gradient regions in terms of second derivatives of certain dependent variables 
may be useful in developing criteria which distinguish the source of 
truncation error from Gibbs' error. Together with the residual data, each 
region having large truncation error can be sorted as to the cause of the 
large truncation error. Criteria for choosing the G weighting in the error 
norms of 3. 2. 3.1 and 3. 2. 3. 4 perhaps can be developed from this basis. See 
Section 5.3 for further discussion of this point. 

At a geometric discontinuity, the sign of the local truncation error 
oscillates at the highest possible frequency of two mesh intervals for an 
ideal difference scheme. This produces a cancellation of the local truncation 
error in the velocity solution. 

It is hypothesized that nonideal difference schemes will exhibit 
two-mesh- interval sign oscillations in the local truncation error estimates 
only at singularities or at .locations which have grid-related problems. 
Otherwise the local truncation error estimates will persist at longer 
wavelengths. It is hypothesized that sums of the same-sign local truncation 
errors are significant to estimating the maximum global error for nonideal 
difference schemes. Useful sums may or may not include the regions of large 
local truncation error depending on the purpose for the error norm. It is 
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hypothesized that the magnitude and the rate change of the local truncation 
error may have use for an error norm where grid juncture in composite grids 
occur or where rapid variations in the dependent variables occur. 

Residual errors and maximum global errors were observed to be directly 
linked. This was examined by computing the discrete continuity balance (local 
mass balance) on each cell. By dividing the local mass balance by the local 
channel cross sectional area, a delta velocity results which, added to the 
local velocity, is the correction necessary to remove the local residual 
error. The maximum global error was reduced to round-off error (below ten to 
the minus ten) when the residual velocity correction was applied successively 
from the entrance region point-by-point through the grid to the exit region. 
Alternatively the maximum global error can be computed directly from the sum 
of the residuals of the same sign divided by the channel cross section at 
which the sign in the residual changes. Control of residual errors is all 
important for satisfying desired global error bounds. It is hypothesized that 
the local residual error should be constrained to some value smaller than the 
local truncation error on the grid which is next to the 'goal grid.' 

5.2 MULTI-GRID 

The form of MG that was used for the computations involves a nonzero 
right-hand side term. With this formulation the discretized continuity 
equation has a mass source right-hand side term which is constructed from the 
estimate of the local truncation error. Fine grid velocity potential data are 
interpolated (restricted) to coarse-grid continuity balances to obtain 
estimates of the local truncation error where global integral is zero for mass 
conservation. Total velocity output that is decoded from solutions of these 
coarse-grid Poisson- type equations are not directly useful (with an academic 
exception). This is a key point about MG output: the total velocity output 

on the coarsest grids may be contaminated with large truncation errors. This 
point is illustrated in Figure 3 for three grid levels. Note that the results 
near the geometric discontinuity are always badly in error. In the coarsest 
grid the local truncation error from the geomtric discontinuity contaminates 
the total velocities at three cell faces where the solution is developed. The 
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extent of the contamination is reduced dramatically as the grid is refined but 

it is only eliminated on the finest grid level where it is exactly zero by 

choice. Any other choice for the finest grid solution would generate worse 

results than that shown; the maximum global error would be larger near the 

discontinuity than occurs in the present example. Therefore the truncation 

(9) 

error extrapolation' cannot be inserted at the finest grid level, only at 
next to the finest grid levels. As shown in Reference 9, it can be used as a 
method for accelerating solution convergence or for generating still finer 
grid solutions (finer than 64 cell cases in the present example) at lower 
cost. Alternatively, a finest grid selection of 32 cells could be used with 
T-extrapolation to get the solution that is shown in Figure 2. 

Standard interpolation fails to be useful for prolongating coarser grid MG 
solutions to finer grid levels. Brandt's FAS-MG formulation is effective for 
this purpose. Estimates of the local truncation error are a direct 
consequence of the FAS-MG process. The maximum value has utility for 
identifying discontinuities. The mean value is a useful guideline of the 
maximum global error. 

It appears desirable to modify conventional applied analysis codes with the 
Brandt FAS scheme so that local truncation error estimates are a routine 
output. This will aid in quickly identifying regions of the analysis domain 
where truncation error problems exist. An optimum MG scheme is not the issue 
for the short term. It is desirable to reduce the labor involved in 
determining where in an analysis domain serious numerical error problems are 
occurring. It may also be feasible to develop error norms that exploit the 
local trunation error estimates of MG so that conventional, semi-adaptive and 
adaptive composite grid technology can achieve high efficiency. 

The grid generation and the PDE solution processes must be drawn together to 
be effective. Composite grid technology should be encouraged. Composite 
grids refer to coupled conformal grids in which nested grids, grid overlays, 
and discontinuous grids are permitted by the analysis approach. 
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Multi “Grid Computer Cost Overhead 

The efficiency of residual error control in terms of computer overhead cost 

for the MG approach was briefly examined. Under the test problem Section 4.3, 

the word "equivalent" sweep was utilized to compare the number of sweeps in an 

MG scheme with FG sweeps. "Equivalent" sweep is defined as follows. 

Neglecting the overhead for the use of the restriction and prolongation 

1 2 

operations for generating RSj^ and estimates, the number of 

iterations on each grid level can be equated to one sweep on the finest grid. 
Thus 16, 8, 4, and 2 are "equivalent" sweeps on grid levels 1, 2, 3, and 4 
respectively relative to level 5 grid. The ratio of the "equivalent" sweeps 
and FG sweeps is not identical with the ratio of (CPU)j^q to (CPU)pj^ To 

achieve parity between these measures of work, a correction factor must be 
empirically generated which corrects the "equivalent" sweeps for the overhead 
of the MG process. This factor is computer and computer code dependent and it 
is not identical with operation counts. No effort has been made to study the 
optimum magnitude of this cost correction factor. For the present 

application, this correction factor is about the same magnitude as the cost to 
sweep the relaxation equation on the finest grid. No optimization of the 
coding was attempted to reduce the size of this factor. Therefore the 
(CPU)[^g/^CP^)fg ■'S approximated by (sweeps) j^g/(sweeps)pg times 

two. A preferred definition of "equivalent" is one that includes this factor. 

Control of the contamination of the total velocity output is correlated with 
the computer work expended in solving the grid equations. The data shows that 
the residual error control efficiency increasingly favors MG over FG as the 
number of grid points is increased. To illustrate this, computations were 
performed as follows. Level 5 grid equations were iterated until a selected 
maximum global error was achieved. The same problem was repeated four more 
times for the maximum global error level using MG. Each problem was 
constrained between two limits of grid level in the MG processes. The grid 
levels are: level 4 - level 5, level 3 - level 5, level 2 - level 5, and 
level 1 - level 5. The (CPU)f^g/(CPU)pQ ratio was estimated by the above 
formula for each problem. It decreased with the extension of the grid level 
separation. The same result was found if the study was performed in the 
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opposite order; namely, (level l)pQ and (level Hevel 2)pg and 

(level 1 - level 2 )(^q, (level 3)pg and (level 1 - level 2 )^^, 

(level 4)pg and (level 1 - level and (level 5)pg and (level 1 - 

level 5)f^Q. This result is in keeping with Brandt's results. For a simple 
elliptic problem, this establishes one type advantage of MG over FG 
procedures: MG is asymptotically more efficient than the FG strategy in 

controlling residual error. Hence the number of grid points that can be 
considered in an analysis with MG is greater than FG for a given computer 
budget. The potential for control of truncation error is thus greater with MG 
than FG strategy for nonideal difference schemes. 

5.3 CONTROL OF THE LENGTH SCALES OF STEEP GRADIENT REGIONS 

The practical utility of the error norms listed in paragraphs 3. 2. 3.1 and 
3. 2. 3. 4 requires quantitative relationships between error norms and such 
parameters as wall skin friction, separation and reattachment points, 
stagnation point location and properties, growth rates of shear layer 
thickness, and the length scale of the resolution of shock waves relative to 
shear layers with which they interact. 

Special attention should be devoted to developing the criteria for the control 

of the length scale of shock waves and stagnation points. One interesting 
(12 13 14) 

recent hypothesis ’ ’ is that shock waves and stagnation points often 

need not be resolved to their true physical length scale. They only need to 

be resolved to length scales of the order-of-magni tude of the key or diffusive 

and boundary regions of the flow field with which they interact. The 

development of quantitative information on the length scale of the large local 

truncation error regions (shock, flame fronts, chemical species fronts, and 

stagnation regions) relative to the diffusive and and boundary regions with 

which they interact is needed. The laminar shock/boundary layer interaction 
(12 13 14) 

problem would be especially suitable for such studies. For 

inviscid computations, singularity length scales must be an input parameter. 

The error norms of equations 3. 2. 3. 1-4, 3. 2. 3. 1-6, 3. 2. 3. 4-1, and 3. 2. 3. 4-3 
have a parameter which must be set to zero in the neighborhood of the 
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singularity regions. The purpose of this parameter is to allow the local 
truncation error to be large relative to the smooth parts of the flow field. 
The detection of these regions is essential. How this can be accomplished is 
a key problem for future work. 

Ideas for detecting these regions can be drawn from mathematical and physical 
features of singularities. For example, shock waves have associated jump 
conditions which characterize the upstream and downstream states in the 
inviscid flow. As the grid refinement process develops, periodic checks can 
be made for the cells that have large second derivatives of these dependent 
variables. Empirical tuning will probably be necessary to define the amount 
of guard mesh around the steep gradient regions that is necessary for limiting 
the grid refinement process. 
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6.0 CONCLUSIONS 


An initial exploratory investigation has been completed toward the development 
of error norms for use as resolution monitors, multi-level adaptive grid 
techniques, and residual error control efficiency for the numerical solution 
of the PDE's of fluid mechanics. Key results are that multi-grid schemes are 
promising as a basis for developing solution resolution monitors, adaptive 
grid techniques, and improved residual error control efficiency. This work 
suggests that multi -grid technology is conceptually straightforward to apply 
to conventional elliptic equation computer codes. Further work is required to 
develop convenient error norms for the local error quantities to guide 
adaptive mesh adjustment with efficient residual error control. 
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7.0 RECOMMENDATIONS 


Various research codes can now be written, based on this initial study, to use 
as vehicles to develop the technology further towards the ultimate goal of 
(semi-) automatic numerical error control in solving the compressible, 
rotational viscous 3-D equations of fluid flow. Reliable numerical error 
monitors for residual and truncation error assessment with efficient control 
methods in conjunction with nested, composite grids should be addressed. The 
following research codes are recommended, with the specific study items listed 
for each of the codes that examine the utility of Ej^, £ 2 , ^Max* ^T» 

^R* ^RT‘ "^Max W* 

3-D Multi grid Potential Code 

A three-dimensional multi-grid code, modeling the full potential equation for 
shockless and shock containing flows, is needed as a test bed for studying 
error norms that relate the maximum global error to the integrals of the local 
multigrid truncation error estimates. With this code, items can be studied 
such as: 

a. Correlation of the maximum global error .with different residual 
tolerances 

b. Convergence criteria for residual tolerance based on maximum local 
truncation error estimates at each grid level or the finest grid 
level 

c. Examine ways to best include steep gradient regions in error norm 
computations. 

1-D Multi -Grid Euler Code 

A one-dimensional multi -grid code, modeling the Euler equations for 
analytically generated channel shapes is needed to study shockless and shock 
containing flows to develop technology features, such as: 

a. The residual control efficiency as affected by the choice of 
relaxation schemes 
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b. The inherent behavior of schemes to remove the acoustic stiffness as 
a stability constraint or modifications needed to achieve the 
infinite speed of sound in low Mach number flow regions 

c. Error norm candidates based on integrals of the local truncation 
error estimate and their relationship to the maximum global error 

d. The inherent capability of multi-grid schemes to remove diffusive 
stiffness from the stability constraints on fine grid levels 

e. Examine approaches to the treatment of steep gradient regions in the 
error norm computations. 

2-D Multi -Grid Euler/Navier- Stokes Code 

A two-dimensional multi-grid code, modeling the Euler and Navier-Stokes 
equations for shockless and shock containing flows is needed to begin applying 
the technology features of error control multi-grid schemes to selected 
problems such as laminar boundary layer and shock/boundary-layer interactions 
on a flat plate. 

Adaptive Embedded Multi -Grid Technology 

The above studies with the various research codes will examine criteria for 
the local multi-grid truncation error estimates that can be utilized for 
labeling cells in which grid nesting is necessary. Recommendations can then 
be made for work tasks to further develop adaptive grid embedding in the 1-D 
Potential, Euler and Navier-Stokes codes and similarly for the 2-D and 3-D 
codes. The design of each of the recommended computer codes should be carried 
out with the goal of flexible grid nesting capability given top priority. 
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CUBIC TRANSITION^ ABRUPT ENLARGEMENT, 31/32 i 2/L 5 1 

SMOOTH CUBIC TRANSITION. 31/32 S Z/L < 1 



Figure 1. One-Dimensional Incompressible Channel Flow 
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Figum 3. Truncation Error Effects on Intermediate Grid Level Solutions 


46 





T/We 



47 








□ 8 Cells on contraction, standard interpolation 

8 Ceils on contraction, Brandt's FAS'InterpoMon 

• 32 Ceils on contraction. Standard fn^pdiatToh 
jQOO 32 Ceils on contraction, Brandt's FAS interpolation. 
Elution exact for 32 ceils on the 
contraction section 
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. Figure 5, Why Brandt's Coarse-to-Fine Grid Correction Equation is Needed 
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